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Abstract. We investigate the slip properties of water confined in graphite-like nano-channels by non¬ 
equilibrium molecular dynamics simulations, with the aim of identifying and analyze separately the influ¬ 
ence of different physical quantities on the slip length. In a system under confinement but connected to a 
reservoir of fluid, the chemical potential is the natural control parameter: we show that two nanochannels 
characterized by the same macroscopic contact angle - but a different microscopic surface potential - do 
not exhibit the same slip length unless the chemical potential of water in the two channels is matched. 
Some methodological issues related to the preparation of samples for the comparative analysis in confined 
geometries are also discussed. 

PACS. 47.ll.-j Computational methods in fluid dynamics, 02.70.-cComputational techniques; simulations 


1 Introduction 

The possibility of slip at the liquid-solid interfaces has 
been debated for over two centuries [T]. The recent blos¬ 
soming of research in micro- and nanofluidics has prompted 
renewed interest in the possibility of slip at the interface 
between a liquid and a solid [51|31|31|SJ[5] , since the classical 
no-slip boundary condition of macroscropic hydrodynamic 
theory does not necessarily hold at those scales [71IH] . In¬ 
stead, the motion of a viscous fluid at a solid interface is 
well described by a partial-slip (Navier) boundary condi¬ 
tion that relates the velocity difference between the solid 
(uw) and the adjacent liquid (ux) to its normal (say along 
z) gradient (see hgure[^ 

( 1 ) 

z=0 

The physical rationale behind the slip phenomenon is a 
balance between the frictional forces and the viscous forces. 
This is apparent in the very same definition of the slip 
length is in a channel of width L (see figure 

is = ry/A - L/2 (2) 

where the viscosity rj = F/ 7 app and the friction coefficient 
A = F/u 

W 0 are determined by the force per unit sur¬ 
face F between the wall and the fluid and the apparent 
shear rate 7app. The role of interfacial slip becomes rele¬ 
vant in conhned micro- and nanofluidic geometries, when 
the surface to volume ratio increases substatially m , so 


. dUj; 

'^W - r\ 

OZ 


Fluid 


Solid 



Solid I 


L 




z 


L 



Fig. 1. A sketch for the steady Couette flow geometry with two 
solid walls moving with opposite velocities FUw A non zero 
slip length (see equation §) can be defined at the boundaries 
when a velocity difference develops between the solid (uw) and 
the adjacent liquid layer (uF)- Consistently, an apparent shear 
rate 7 app F 2uwlL develops. 


that the slip length may become comparable with the sys¬ 
tem size, is ~ L: an accurate understanding of nanoscale 
friction phenomena at the fluid-solid interfaces is indeed 
paramount to the design of micro- and nanofluidic de¬ 
vices aimed at optimizing mass transport against an over¬ 
whelming dissipation barrier [ini[iii0. In the recent years, 
progress has been made in providing a coherent descrip¬ 
tion of slip motion close to solid boundaries, both from 
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the experimental m and theoretical [l2] side, and gen¬ 
eral consensus on the occurrence of slip on partially wet¬ 
ting substrates has emerged. Nevertheless, the physical 
mechanisms leading to fluid slip are still not completely 
understood [HlHlIillll]. In a recent paper, Huang and 
coworkers m, using molecular dynamics simulations of 
an atomistic water model, studied the interfacial hydro- 
dynamic slippage at various hydrophobic surfaces. The 
measured slip lengths range from nanometers to tens of 
nanometers and were shown to almost collapse onto a 
single curve as a function of the static contact angle 9 
characterizing the surface wettability m, thereby sug¬ 
gesting a quasi-universal relationship ig oc (1 -I- cos0)~^ 
[H], with 9 the equilibrium contact angle formed by liq¬ 
uid droplets at the solid surface. Rationalizing this picture 
in terms of the statistical properties of atomistic motion 
close to solid boundaries is a challenging task [T^nfT^ITni 
EolEI], especially when different molecular mechanisms of 
slip may coexist depending on the external driving 

force (the applied shear), liquid atoms may hop from one 
equilibrium site to another of the liquid-solid energy land¬ 
scape, or even display a collective motion of entire layers 
of atoms slipping together. More recent numerical simu¬ 
lations m even suggested that the dichotomy between 
hydrophobicity and large slip might be purely coinciden¬ 
tal and that hydrophilic surfaces can show features typ¬ 
ically associated with hydrophobicity. The description of 
the fluid and surface at atomistic level is essential for these 
investigations, as small changes in the surface properties 
can lead to surprisingly different structural and dynam¬ 
ical properties of the fluid, as it is seen for example for 
hydroxylated surfaces Haul. 

Another issue is that of the shear-rate dependence of 
the slip length: Thomson and Troian discovered, for exam¬ 
ple, that the slip length can follow a universal curve as a 
function of the applied shear rate, leading to a divergence. [25] 
Priezjev, however, showed that an important role in this 
situation is played by static surface roughness |26|. As it 
turned out, displacements of the atoms in the surface layer 
as small as those generated by thermal fluctuations in 
graphite (about one tenth of an Angstrom) are enough 
to remove the slip length divergence, demonstrating how 
sensitive this quantity is on the microscopic detail of the 
surface [27]. 

Further analysis on the molecular mechanisms respon¬ 
sible for these stimulating results is needed: the picture 
of a quasi-universal relationship between wettability and 
slip length is surely appealing, but it should be born in 
mind that as a dynamical quantity, the slip length can 
not depend only on static properties like the contact angle 
(that is, on the interfacial free energies). Indeed, a quan¬ 
tity with the dimension of a diffusion coefficient appears 
in the approximate expression for the slip length obtained 
from Green-Kubo relations [I^l^ . More detailed calcula¬ 
tions suggest this quantity to be related to the collective 
solvent diffusion tensor (29). This clarifies the meaning of 
the quasi-universal relation proposed in m , which strictly 
speaking holds true only if the dynamical properties of the 
solvent in interaction with the different substrates do not 
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Fig. 2. Top: side view of a detail of the three atomic layers 
composing the random quenched surface. Atoms with purely 
repulsive Lennard-Jones interaction are darker (red), while 
those with the enhanced attractive part of the potential are 
lighter (gray). On top of the atomic layers the equipotential 
surfaces at energy E = kT and E = 2kT are presented (not 
to scale in the direction perpendicular to the surface), to show 
the shape of the attractive basin next to the surface. Bottom: 
the top view of the detail of the random quenched surface. The 
atoms are depicted as in the top panel, and a translucent rep¬ 
resentation of the equipotential surface E = kT is overlaid, in 
order to show the potential energy inhomogeneity. 


change much depending on the substrate nature. The pres¬ 
ence of marked surface inhomogeneities, the possibility of 
making hydrogen bonds, or the presence of dissociated 
charged groups at the solid-liquid interface can of course 
drastically change this picture. 

In this paper we address this issue by extending the 
systematic approach that allowed us to tackle the slip 
length divergence problem [2 7j. The approach consists in 
employing model surfaces with different microscopic prop¬ 
erties (i.e. a different interaction potential between solid 
wall particles and water), but the same macroscopic con¬ 
tact angle of the solvent. In our previous work, however, 
we used the same number of water molecules in each of 
the simulations with different surfaces, hence, water was 
simulated at different chemical potentials. Here, we inves¬ 
tigate the effect of this difference in chemical potential 
on the slip length, and show that two different surfaces, 
even if they share the same macroscopic contact angle, can 
indeed lead to a considerably different slip lengths. This 
difference, however, vanishes when the chemical potential 
of water in the two channels is matched. This has some 
important implications from the practical and from the 
theoretical point of view, which are here discussed. 


2 Methods and Systems 

Using an in-house modified version of the Gromacs suite |30j. 
we simulated a slab of water molecules, modelled with the 
SPG/E potential [3T] and confined between two graphite¬ 
like walls, each consisting of three atomic layers. The gap 
between the water-exposed layers is 6.982 nm, and the box 
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edge sizes are 13.035, 16.188 and 8.0 nm in the x, y and 
z directions, respectively. 

The electrostatic interaction was calculated using the 
smooth particle mesh Ewald method [32], together with 
Yeh-Berkowitz correction [33] to remove the contribution 
from periodic images along the wall normal. Time integra¬ 
tion of the equations of motion was performed using the 
leapfrog algorithm |34] with an integration time step of 1 
fs. Water molecules are kept rigid using the SETTLE algo¬ 
rithm |35| . Both in equilibrium and non-equilibrium simu¬ 
lations, a Nose-Hoover thermostat, coupled to the degrees 
of freedom not subjected to external forces, was employed 
to keep temperature constant. The atoms in the substrate 
are kept fixed in the equilibrium MD simulations, and are 
moved at a fixed velocity in the NEMD simulations. No 
center of mass removal procedure was applied. 

We investigated two different setups, with respect to 
the wall-liquid interactions. In a first setup (here identi¬ 
fied as “standard”) the interaction potential of wall atoms 
with oxygen atoms is given by a Lennard-Jones poten¬ 
tial C7(r) = C'i 2 r“^^ — Cgr”®, with Cq = 2.47512 x 10“^ 
kJ/nm® and C 12 = 2.836328 x 10“°® kJ/nm^^|5Bj up to 
a distance of 0.9 nm. Above that distance, the potential 
was smoothly switched to zero at 1.2 nm, using 4th or¬ 
der polynomials. The second setup consists of a “random 
quenched” functionalization of the previous one, realized 
by making 40% of wall atoms purely repulsive {C% = 0), 
and increasing the interaction strength of the remaining 
ones by a factor a, to be tuned in order to achieve the 
same macroscopic contact angle as in the standard case. 
Equipotential energy surfaces and a snapshot of the three 
layers are shown in Figj^for the random quenched case. 

The procedure we used to investigate the slip length 
for different surfaces removing the biasing dependence on 
both contact angle and chemical potential, can be summa¬ 
rized as follows: a) we first generated the depleted surface 
by removing the attractive term of the Lennard-Jones in¬ 
teraction for a random selection of surface atoms; b) we 
calibrated the Lennard-Jones interaction strength a of the 
remaining atoms to obtain the same macroscopic contact 
angle of the standard surface case, using the generalized 
Young equation to extrapolate the contact angle to its 
macroscopic limit for a sessile droplet; c) we determined 
the number of water molecules that yields equal chemical 
potential for the systems with standard and depleted sur¬ 
faces in a slit-pore configuration; d) we calculated the ve¬ 
locity profile induced by moving at constant speed and in 
opposite direction the two slabs of the slit pore, at a shear 
rate low enough to be in the constant slip length regime ; 
e) we eventually extracted the apparent slip length from 
a linear fit of the velocity profile, performed in the central 
part of the slit pore. 

3 Results and discussion 

The contact angle. Given the known dependence of the 
slip length on the contact angle [T3] , an equal wetting is a 
prerequisite to compare the slip length of water on two mi¬ 
croscopically different surfaces and understand how much 


the slip length is affected by different types of surface in¬ 
homogeneities. The actual value of the factor a has been 
chosen so that water wets the two surfaces with a compa¬ 
rable macroscopic contact angle. To compute the macro¬ 
scopic contact angle, we employed the procedure suggested 
by Werder and colleagues |37|. For each surface type, three 
droplets of different sizes (about 1.8 x 10®,14 x 10® and 
34 X 10® water molecules, respectively) on six monoatomic 
layers were simulated starting from a rectangular droplet 
shape, which is equilibrated for 500 ps (shape relaxation 
is observed to take place within few tens of ps). Density 
profiles along the radial direction within slabs at different 
heights z, were sampled from the next 1 ns of simulation. 

A best fit to a sigmoid function is then employed to iden¬ 
tify the location of the Gibbs dividing surface Rg{z). A 
best fit to a circumference is eventually performed using 
the points of the locus Rg{z), to extract the droplet radius 
i?, base radius Rb and contact angle 9. The extrapolation 
of the contact angle to infinite is then performed via 
the linear fit o the generalized Young equation 

cos( 6 i) = cos( 6 'oo)-— (3) 

iLv Rb 

where r is the line tension, cos(doo) = isv — ISL, and 
7 iy, 75 y and 75 are the liquid-vapour, solid-vapour and 
solid-liquid surface tensions, respectively. The lateral size 
of the substrate layers is chosen to be at least a factor 1.5 
larger than the droplet base radius. 

The calibration of the Lennard-Jones interactions turned 
out to be a delicate procedure. The cosine of the macro¬ 
scopic contact angle, reported in Figj^ presents a rather 
pronounced scattering, arising from the fitting of the ra¬ 
dius Rb and of the angle 0 , as well from the extrapolation 
to macrosopic droplet sizes, which limited the accuracy 
in the determination of the macroscopic contact angle to 
about ±5 deg. This limitation, however, turned out not 
to be overly restrictive for our purposes, as the contact 
angle for the standard surface is about 40 deg, and the 
slip length should show only a weak dependence on the 
contact angle, for low values of the latter. In this case, the 
deviation of the slip length due to a contact angle change 
of 5 deg would be roughly 6 %. As an outcome of this cali¬ 
bration procedure, we have chosen the interaction of that 
60% of atoms that still retain the attractive part of the 
potential to be a factor a = 1.977 stronger than in the 
case of the standard surface. 

The cosine of the macroscopic contact angle, as it is 
seen in Figj^ shows a linear dependence as a function of 
the scaling factor applied to the Lennard-Jones interaction 
of the attractive sites. A mean field approach is enough 
to explain this behaviour, as it is known that the liquid- 
solid surface tension in a fluid that interacts through the 
Lennard-Jones potential with the surface depends linearly 
on the Lennard-Jones energy]!^. It is worth noticing that 
also the values of the line tension r, although characterized 
by a high spread, can still be described reasonably well by 
a linear dependence (reduced = 0 . 8 ) on the interaction 
strength. The line tensions ranged from negative values 
(r ~ -30 kJ/mol/nm for the smallest interaction strength), 
to positive ones (r ~ 80 kJ/mol/nm). 
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Fig. 3. Top: radial profiles of three different water droplets on 
the standard surface. The large fluctuations next to the top 
of the droplet are caused by the small volume element, and 
consequent low number of molecules associated to it. Bottom: 
simulation snapshot of one of the droplets, showing also some 
molecules in the vapour phase. 


The chemical potential. Before being able to carry out 
the Couette flow simulation for the modified surface, the 
proper water content between the two slabs had to be 
determined. The natural thermodynamic control parame¬ 
ters for a fluid confined between two flat interfaces are the 
accessible volume V, the temperature T, and the excess 
chemical potential of the fluid, since the latter is usu¬ 
ally at thermal and chemical equilibrium with a reservoir. 
A series of MD simulations of a slit pore with different wa¬ 
ter content were hence performed, in order to identify the 
number of water molecules JV^ that leads to the same ex¬ 
cess chemical potential of the standard surface case. The 
chemical potential was computed using the Widom inser¬ 
tion method |39| in its variant for inhomogeneous systems 
m- The chemical potential of a system with water 
molecules at one point r in space is given by 


1, |’(exp(-MC/(r))> 

= —Se- 


( 4 ) 



LJ energy ratio 

Fig. 4. Cosine of the macroscopic contact angle as a func¬ 
tion of the ratio between the Lennard-Jones parameters e = 
(71/(4(712) of the atoms from set S 2 and that of the original 
force-field parameters. The value of the macroscopic contact 
angle for the standard surface is marked as a horizontal line. 
The circle on the point at a ratio of 1.977 marks the case chosen 
for the subsequent steps. 



Fig. 5. Chemical potential of water conhned between two slabs 
(circles, standard substrate; squares, depleted substrate). Solid 
lines are the result of a best fit to a quadratic function. 


where AU is the change in potential energy due to the 
insertion of an additional molecule at position r with ran¬ 
dom orientation, and the canonical average was performed 
by sampling equilibrium configurations of the system with 
Nn, water molecules. For the Widom method, water molecules 
were inserted in those region where the density of water 
ensured proper statistics, namely, not too close to the sur¬ 
face. At chemical equilibrium, ^Xx{y) has to be equal at ev¬ 
ery point in space, and this can be exploited to sample the 
chemical potential with high accuracy and relatively lit¬ 
tle computational effort. In the slab system, translational 
invariance along the x and y directions can be assumed 
- not too close to the surface - even if the surface is not 
homogeneous, so that the test molecule can be inserted at 
random positions on planes at different heights z. For the 
chemical potential calculation, the systems with different 
water content were simulated at equilibrium for 40 ns, sav- 
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Fig. 6. Slip length of water, as a function of the imposed shear 
rate, for three different systems. Circles: standard surface with 
water density pi. Squares: random quenched surface with water 
density pi. Triangles: standard surface with water density p 2 
(and with chemical potential matching that of squares). 

ing configurations every 1 ps for the Widom test particle 
insertion analysis. For each frame stored, 10"^ positions in 
a plane parallel to the substrate surface were randomly 
chosen as insertion sites for a test water molecule, and the 
energy difference was sampled according to the Widom 
scheme to compute the chemical potential. The obtained 
chemical potentials are reported in Figj^ together with 
the result of a fit to a quadratic function. The chemical 
potentials of the two systems, for a water content lower 
than 44 X10^ differ by about 2 kJ/mol, but become compa¬ 
rable in the proximity of — 45 x 10^. Since in this work 
we did not use any third reference system (open reservoir), 
the water content of one system can be chosen arbitrarily, 
and that of the other system can be changed to match the 
chemical potential. We hence decided to use = 43413 
for the standard surface, and, making use of the quadratic 
interpolation, we identified the value that yields the same 
chemical potential = —11.4±0.12 kJ/mol) for the de¬ 
pleted surface, namely, = 44133. The two correspond¬ 
ing densities will be denoted as pi and p 2 , respectively 
(P2 > Pi)- 


The slip length. Slip lengths are eventually calculated 
by extrapolating the velocity profile of a Couette flow in¬ 
duced by imposing constant velocity on two parallel, iden¬ 
tical slabs made of the atoms of the standard or depleted 
surface. Configurations are saved every 1 ps to analyze 
the velocity profile {x component) along the z direction 
and to extract the apparent slip length, using a linear 
fit of the resulting Couette flow. In Figj^we show the slip 
length obtained for three different systems at various shear 
rates. The three systems are: (i) the random quenched 
surface channel with water density pi (squares); (ii) the 
standard, unmodified surface with either (ii) pi (circles) 
or (iii) p 2 (triangles) For all three cases, the slip length 
starts increasing noticeably at shear rates higher than 


Fig. 7. Slip length of water in the pore with standard sur¬ 
face, as a function of the chemical potential (from the best fit), 
obtained for an imposed shear rate of 0.0047 ps“^. The water 
content N-w of the channel is reported on the plot next to each 
point. The lowest and highest water content correspond to the 
densities pi and p 2 , respectively. At the density p 2 the chem¬ 
ical potential of the standard surface system matches that of 
the random quenched surface system with water density pi 
(square) 

~ 0.01 ps“^, but is shear-independent below shear rates of 
about 0.005ps“^. In this plateau region the systems with 
same water content but different surface potential show a 
noticeably different slip length. On the contrary, match¬ 
ing the chemical potential drives the system to the same 
slip length (always in the plateau region). This result, we 
would like to stress, is valid under the condition that the 
macroscopic contact angle for water droplets is the same 
for both surfaces. 

By calculating the slip length for the systems with the 
standard surface at different water content it is possible 
to investigate quantitatively the dependence on the slip 
length on the excess chemical potential, as shown in Fig.[^ 
Even though the statistical uncertainty is rather large, one 
can observe a definite trend. The magnitude of the slip 
length change is particularly important, as it passes from 
9 to 12.5 nm (an increase of roughly 30%, compared to 
a decrease of about 15% in chemical potential). To pro¬ 
vide some practical terms of comparison, one can consider 
that the temperature coefficient for the chemical poten¬ 
tial of water at standard conditions is dp/dT = —69.9 
,I/mol/K|41). hence, an equivalent change of 2 kJ/mol in 
the chemical potential of water could be realized, for ex¬ 
ample, by raising its temperature by about 30 K. 

It is an important question, wheter a change from 9 to 
12.5 nm is a significant one. In a Poiseuille (cylindrical) 
flow setup, the relative rate of work WijW^ necessary to 
achieve a target flux Q in two systems characterized by 
the slip lengths l\ and respectively, is 

W, f l+4£2/R Y 

VFa [i+4£i/eJ ’ ^ ^ 

where R is the tube radius. This result can be obtained 
considering that the work rate (per unit length) in pres- 
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ence of a slip velocity u{R) is 


w = 


1 - 


ttR^u{R) 

Q 


( 6 ) 


where dp/dz is the pressure gradient. The first term on the 
right hand side comes from the fluid deformation, while 
the second one from the friction at the boundary [35], and 
can be cast in this form noticing that the friction force at 
the surface has to balance the force generated by the pres¬ 
sure gradient, 2ttR< 7 = TrR^dp/dz, where a is the stress 
tensor. For a Poiseuille flow, the slip velocity is 

= (7) 

and the relation between pressure difference and flux is 




dp ttR^ 
dz 8r] 



( 8 ) 


from which, after a bit of algebra, Eq. can be derived. 
With the present values of slip length, £i = 12.5, = 9, 

and using R = L, the outcome is a signihcant reduction, 
Wi ~ 0.54bF2, of the work rate necessary to sustain the 
flow. 
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